Pricing and Hedging Asian Basket Options 
with Quasi-Monte Carlo Simulations 



Nicola Cufaro Petroni 
Dipartimento di Matematica and TIRES, Universita di Bari 
INFN Sezione di Bari 
via E. Orabona 4, 70125 Bari, Italy 
email: cufaro@ba.infn.it 

Piergiacomo Sabino 
Dipartimento di Matematica, Universita di Bari 
via E. Orabona 4, 70125 Bari, Italy 
email: sabino@dm.uniba.it 

July 17, 2009 
Abstract 

In this article we consider the problem of pricing and hedging high-dimensional 
Asian basket options by Quasi-Monte Carlo simulation. We assume a Black-Scholes 
market with time-dependent volatilities and show how to compute the deltas by 
the aid of the Malliavin Calculus, extending the procedure employed by Montero 
and Kohatsu-Higa [lj. Efficient path-generation algorithms, such as Linear Trans- 
formation and Principal Component Analysis, exhibits a high computational cost 
in a market with time-dependent volatilities. We present a new and fast Cholesky 
algorithm for block matrices that makes the Linear Transformation even more 
convenient. Moreover, we propose a new-path generation technique based on a 
Kronecker Product Approximation. This construction returns the same accuracy 
of the Linear Transformation used for the computation of the deltas and the prices 
in the case of correlated asset returns while requiring a lower computational time. 
All these techniques can be easily employed for stochastic volatility models based 
on the mixture of multi-dimensional dynamics introduced by Brigo et al. [21 13] ■ 

1 Introduction and Motivation 

In a recent paper Dahl and Benth j4] have investigated the efficiency and the 
computational cost of the Principal Component Analysis (PCA) used in the Quasi- 
Monte Carlo (QMC) simulations for the pricing of high-dimensional Asian basket 
options in a multi-dimensional Black-Scholes (BS) model with constant volatilities. 
In particular they have shown the essential role of the Kronecker product for a 
fast implementation as well as for the analysis of variance (ANOVA) in order to 
identify the effective dimension (see below). Indeed the convergence rate of the 
QMC method is O (^N^ 1 log d N^j , where N is the number of simulation trials, and 
d the nominal dimension of the problem. This implies that the theoretically higher 



1 



N Cufaro Petroni and P Sabino: Asian Basket Options and QMC 



2 



asymptotic convergence rate of QMC could not be achieved for practical purposes 
in high dimensions. On the other hand, some applications in finance (see Paskov 
and Traub [5]) have shown that QMC provides a higher accuracy than standard 
Monte Carlo (MC), even for high dimensions. 

To explain the success of QMC in high dimensions Caflisch et al. [5] have 
introduced two notions of effective dimensions based on the ANOVA of the inte- 
grand function. Consider an integrand function / and a MC problem with nominal 
dimension d. Let A — {1, . . . , d} denote the labels of the input variables of the 
function /: then the effective dimension of /, in the superposition sense, is the 
smallest integer ds such that J2\ u \<d s (j2 (/«) — P°' 2 (/)' where f u is a function 
with variables in the set u C A, c 2 (-) denotes the variance of the given function, 
|u| is the cardinality of the set and < p < 1, for instance (p = 0.99). The effec- 
tive dimension of /, in the truncation sense is the smallest integer dr such that 
SuC{i 2 d T } a2 (S u ) = P a2 (f)- Essentially, the truncation dimension indicates 
the number of important variables which predominantly capture the given func- 
tion /. The superposition dimension takes into account that for some integrands, 
the inputs might influence the outcome through their joint action within small 
groups. 

The PCA decomposition only permits a dimension reduction without taking 
into account the particular payoff function of a European option. In contrast, Imai 
and Tan [7] have proposed a general dimension reduction construction, named Lin- 
ear Transformation (LT), that depends on the payoff function and that minimizes 
the effective dimension in the truncation sense. They have shown that the LT ap- 
proach is more accurate than the standard PCA, but has a higher computational 
cost. In a previous paper one of the authors [8] has discussed how to implement 
this technique quickly and - with a slower computer - has obtained computational 
times that are about 30 times smaller that those originally presented by Imai and 
Tan 0. 

In the present study we consider time-dependent volatilities: as a consequence it 
is not possible to rely on the properties of the Kronecker product, and the problem 
is computationally more complex. In order to simplify the computational complex- 
ity, we present a fast Cholesky (CH) decomposition algorithm tailored for block 
matrices that remarkably reduces the computational cost. Moreover, we present 
a new path-generation technique based on the Kronecker Product Approximation 
(KPA) of the correlation matrix of the multi-dimensional Brownian path that re- 
turns a suboptimal ANOVA decomposition with a substantial advantage from the 
computational point of view. 

Our numerical experiment consists in calculating the Randomized QMC (RQMC) 
estimation of the prices and the deltas of high-dimensional Asian basket options 
in a BS market with time-dependent volatilities. In order to compute the deltas, 
we extend to a multi-assets dependence the procedure employed by Montero and 
Kohatsu-Higa [T] in a single asset setting. This procedure is based on the Malliavin 
Calculus and allows a certain flexibility that can enhance the localization technique 
introduced by Fournie et al. [9] . As far as the computation of Asian options prices 
is concerned, the KPA and LT approaches are tested both in terms of accuracy 
and computational cost. We demonstrate that the LT construction becomes more 
efficient than the PCA even from the computational point of view, provided we 
use the CH algorithm that we present and the approach described in [8] . The KPA 
and the PCA constructions perform equally in terms of accuracy, with the former 
one requiring a considerably shorter computational time. However, the KPA and 
the LT display the same accuracies in the computation of the deltas. Moreover, we 
compare our simulation experiment - also using the standard CH and the PCA de- 
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composition methods - with pseudo-random and Latin Hypercube Sampling (LHS) 
generators. 

Remark finally that all the methods described here can accommodate a market 
with stochastic volatility where the evolution of the risky securities is modeled by 
a mixture of multi-dimensional dynamics as in the papers by Brigo et al. [H [3] • It 
is noteworthy to say that none of these constructions can be applied to Heston-likc 
multi-dimensional stochastic volatility models. In principle, we might still use the 
LT for the Euler discretization of the Heston model, but this could be no longer 
applicable with more realistic schemes that involve discrete random variables as 
proposed for instance by Alfonsi [TU] . 

The paper is organized as follows: Section [2] describes Asian options. Section 
[3] discusses some path-generation techniques and in particular, presents the fast 
CH algorithm and the KPA construction. Section |4] shows the numerical tests for 
the Asian option pricing. Section [5] explains how to represent the deltas of Asian 
basket options as expected values with the aid of Malliavin Calculus and shows the 
estimated values by RQMC. Section [5] summarizes the most important results and 
concludes the paper. 



2 Asian Basket Options 

Assume a multi-dimensional BS market with M risky securities and one risk-free 
asset. Denote B (t) = (Bi (t) , . . . , Bm (t)) an M-dimensional Brownian motion 
(BM) with correlated components and (!Ft)t>o the filtration generated by this BM. 
Moreover, denote pik the constant instantaneous correlation between Bi(t) and 
■Bfc(t), Si (t) the i-th asset price at time t, Oi (t) the instantaneous time-dependent 
volatility of the i-th asset return and r the continuously compounded risk-free rate. 
In the risk-neutral probability, we assume that the dynamics of the risky assets are 



dS l (i) = rSi (t) dt + en (t) S, (t) dBi (t) 
The solution of Equation ([I]) is 



1. 



,Af. 



(1) 



Si{t) = Si(0)exp 



ds 



Oi (s) dBi (s) 



i = l,...,M. (2) 



Discretely monitored Asian basket options are derivative contracts that depend on 
the arithmetic mean of the prices assumed by a linear combination of the underlying 
securities at precise times t\ < t<i • ■ • < — T, where T is the maturity of the 
contract. By the risk-neutral pricing formula (see for instance Lamberton and 
Lapeyre [11]) the fair price at time t of the contract is 



a (t) = e r(T -^E 



M N 

J2^2wij S l (t J ) - K 

= 1 3=1 



(3) 



with the assumption that ^ ■ Wij = 1. 

Pricing Asian options by simulation hence requires the discrete averaging of the 
solution @ at a finite set of times {t\ 1 . . . ,t^}- This sampling procedure yields 



S^) = S' l (0)exp 



' r - ) dt + Zi(t jt 



l,...,M,j = l,...,N, 



(4) 
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where the components of the vector 
(Zi(ti),...,Zi(tjv); Z 2 (ti),.. 



> Z 2 (ijv); • • • ; ZM(h), . . . , Zm(£aO) J 



are M x N normal random variables with zero mean and the following covariance 
matrix 

/ £(ii) S(t!) ... E(ti) \ 
E(ti) £(t 2 ) ... £(i 2 ) 



Emjv = 



V E(ti) S(t 2 ) 



(5) 



S(tjv) / 



where the elements of the MxM submatrices £(i„) are (E(t„)) ifc = J™ pik<Ti{s)<7k{s)ds 
with i,k = 1, . . . , M; n = 1, . . . , N. This setting would allow time-dependent cor- 
relations as well. In the case of constant volatilities the covariance matrix is 



Emw — 



( tiE tiE ... tiE ^ 



V *iS t 2 £ 



(6) 



t^E y 



where now E denotes the M x M covariance matrix of the logarithmic returns of 
the assets. It follows from the last equation that the covariance matrix Y,mn can 
be represented as R <g> E, where ® denotes the Kronecker product and R is the 
auto-covariance matrix of a single BM. This simplification is not possible in the 
case of time-dependent volatilities. We recall that the elements of R are 

Ri„=tiAt n , l,n=l,...N. (7) 

R has the peculiarity to be invariant for a reflection about the diagonal. 

Definition 1 (Boomerang Matrix). Let B 6 R"b x "b be a square matrix and let 
b = (bi , . . . , b nB ) G R nB . B is a boomerang matrix if 



Bh P = bhA P , h,p = 1, . . . ,hb- 
We call b the elementary vector associated to B. 



(8) 



As a consequence R is boomerang, and in general the auto-covariance matrix of 
every Gaussian process is boomerang. This definition can also be extended to block 
matrices as follows. 

Definition 2 (Block Boomerang Matrix). Partition the rows and the columns of 
a square matrix B £ R™b x,1 b to obtain: 



( B, 



B = 



Sip \ 



(9) 



y Bpi . . . Bpp ) 

DxD designates the (h,p) square submatrix and 



where for h,p — 1, . . . , P, Bh p G 
n B = P x D. Given P matrices B\, . . . , B P with B h G R DxD , h 
a boomerang block matrix if 



1.....P, B is 



B, 



hp 



B h 



h,p = l,...,n B - 



(10) 



We call b = (B\, . . . ,Bp) T the elementary block vector associated to B. 



N Cufaro Petroni and P Sabino: Asian Basket Options and QMC 



5 



From these definitions we have that Ea/tv is block boomerang. 

The payoff at maturity of the Asian basket option now is a(T) = (,g(Z) — K) + 
with 

MxN 

g(Z) = exp( Mfc + Z fe ) (11) 
fc=i 

where Z ~ Af(0, Sjwjv) arid 

f tk2 <?l (*) 

H k = \n{w klk2 S kl (0)) + rt k2 - / (ft (12) 

with fci = (k - 1) mod M; fc 2 = L( fc - 1) /MJ + 1; = 1, . . . , M, where |xj denotes 
the greatest integer less than or equal to x. 



3 Path-generation Techniques 

From the previous discussion it comes out that the pricing of Asian basket options 
by simulation requires an averaging on the sample trajectories of an M-dimensional 
BM. In general, if Y ~ A/"(0, Sy) and X ~ A/"(0, 1) are two iV-dimensional Gaussian 
random vectors, we will alawys be able to write Y = CX, where C is a matrix 
such that: 

Sy = CC T . (13) 

and the core problem consists in finding the matrix C . In our case Y,y coincides 
with Emtv of Equation ([5]). The accuracy of the standard MC method does not 
depend on the choice of the matrix C because the order of the random variables 
is not important. However, a choice of C that reduces the nominal dimension 
would improve the efficiency of the (R)QMC. In the following we discuss some 
possibilities. 



3.1 Cholesky Construction 



The CH decomposition simply finds the matrix C among all the lower triangular 
matrices. In the case of constant volatilities the matrix Emtv is the Kronecker 
product of R and S, and the Kronecker product shows compatibility with the CH 
decomposition (see Pitsianis and Van Loan [E]). Indeed, denote Cs MN , Cr and 
Cs the CH matrices associated to Emat, R and S respectively; we then have 



(14) 



This now allows a remarkable reduction of the computational cost: it turns out in 
fact that a O ((M x N) 3 ) computation is reduced to a O (M 3 ) + O (N 3 ) one. 

When time-dependent volatilities are considered, however, we can no longer 
use these properties of the Kronecker product. In any case, since Smjv is a block 
boomerang matrix, we can use the following result: 



Proposition 1. Let B £ 

where B) t 6 
block vector. Cb 



j xn B a block boomerang matrix and let (B\, . . . , Bp) 
DxD f h — 1, . . . ,P with ns = P x D, be its associated elementary 
the CH matrix associated to B, is given by: 



C B = 



( Ci 
: C 2 



\ 



(15) 



V C 1 C 2 



c P J 
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where the D x D blocks Cu, h = 1, . . . , P are 

C h = Ch.ol(B h -B h _ 1 ) (16) 

with Choi denoting the CH factorization, and we assume Bo = 0. 

Proof. Consider the h th row of Cb and the m th row of its transposed matrix; we 
then have 



(c 1 ,...,c^o,...,o) T - (cf[,...,cl,o,...,o) = ^qcf 

hf\m 

= ^(Bi - Bi-i) = B hAm 
and this concludes the proof. □ 
3.2 Principal Component Analysis 

Acworth et al. [13] have proposed a path generation technique based on the PCA. 
Following this approach we consider the spectral decomposition of Smtv 

S M jv = EAE T = (£A 1/2 )(£A 1/2 ) T , (17) 

where A is the diagonal matrix of all the positive eigenvalues of Emjv sorted in 
decreasing order and E is the orthogonal matrix (EE T = I) of all the associated 
eigenvectors. The matrix C solving Equation (fT3)) is then EK 1 / 2 . The amount of 

variance explained by the first k principal components is the ratio: £^ =1 * where 

d is the rank of Smjv- The PCA construction permits the statistical ranking of the 
normal factors, while this is not possible by the CH decomposition. For the market 
with constant volatilities, the Kronecker product reduces this calculation into the 
computation of the eigenvalues and vectors of the two smaller matrices R and 
S. All these simplifications are no longer valid for the time-dependent volatilities. 
However, we can reduce the computational cost for the PCA decomposition in the 
following way. 

Take Mi, Ma, M3 and M4, respectively p x p,p x q,q x p and q x q matrices, 
and suppose that Mi and M4 are invertible. Assume 



M = 



Mi M 2 
M 3 M 4 



and define Si = M 4 - MsM^Ma and S4 = Mi - MjM^M?, the Schur com- 
plements of Mi and M4, respectively. Then by Schur's lemma the inverse M _1 

M -i_( S, -M^M 2 S^\ 

Taking into account the previous result it is possible to prove the following propo- 
sition 

Proposition 2. Let B £ R" bX " b be a block boomerang matrix and let (B\, . . . , Bp) T , 
where Bh G M £,xI? , h = 1, . . . , P with ns = P x D , be its associated elementary 
block vector. The inverse of B is symmetric block tri-diagonal. The blocks on 
the lower (and upper) diagonal are 7j = — (£>f+i — B{) 1 , I = 1, . . . , P — 1 while 
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those on the diagonal are D m — (B m — _B m _i) 1 (B m+ i — B n 
m = 1, . . . , P, with the assumption that B — -B/v+i = 0: 



0(B 



m+1 



B m ) 



( D 1 
Ti 


V o 



D 2 
T 2 







2>_! 



\ 





Tp-x 
Dp 



(19) 



This property can be used to reduce the computational cost of evaluating the 
PCA decomposition in the case of time-dependent volatilities and in general for 
multi-dimensional Gaussian processes. Indeed, if B is a non-singular square matrix 
then the eigenvalues of the B~ l are the reciprocal of the eigenvalues of B and the 
eigenvectors coincide. 



3.3 Linear Transformation 

Imai and Tan [7] have considered the following class of LT as a solution of Q13p : 

c lt = c ch A ( 20 ) 

where C ch is the CH matrix associated to the covariance matrix of the normal 
random vector to be generated, and A is an orthogonal matrix, i.e. AA T = I. 
The matrix A is introduced with the main purpose of minimizing the effective 
dimension of a simulation problem in the truncation sense. Imai and Tan [7] have 
proposed to approximate an arbitrary function g, such that (g — K) + is the payoff 
function of a European derivative contract, with its first order Taylor expansion 
around e 

n n 
\ - ag 



= 9(e) + 



^ dei 

i=i 



(21) 



The approximated function is linear in the standard normal random vector Ae. 
Considering an arbitrary point about which we form the expansion, such as e = 0, 
we can derive the first column of the optimal orthogonal matrix A* . It is possible to 
find the complete matrix by expanding g about different points and then compute 
the optimization algorithm. Imai and Tan [7] have set: t\ = = (0, 0, ... , 0), e 2 — 
(1, 0, . . . , 0), . . . , e„ = (1, . . . , 1, 0), where the fc-th point has k — 1 non-zero compo- 
nents. The optimization can then be formulated as follows: 



max 
a k eR n 



dg 



k = 1, 



.,n, 



(22) 



subject to ||A.k|| = 1 and A*j • A.k = 0; j = 1, . . . , k — 1; k < n. In the case of 
Asian basket options we have 



NAI 



9(e) = 9(e) + J2 



1=1 



NM / NM \ 

^ exp I jjbj + C ik e k C a Ae ; . 

i=l \ k=l ) 



(23) 



Imai and Tan [JJ have proved the following result: 

Proposition 3. Consider an Asian basket options in a BS model, define: 



d ( P ) 
B (p) 



3 (w/N + ELl C MN,k) 



(C Ch ) T (d(p)), p = l,...,MN. 



(24) 
(25) 
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Then the p-th column of the optimal matrix A* is 

B (p) 
]BW|| 



A *p = ±^57^T P = 1,---,MN. (26) 



The matrices C* k , k < p have been already found at the p — 1 previous steps. A. p 
must be orthogonal to all the other columns. This feature can be easily obtained 
by an incremental QR decomposition as described in Sabino [5]. 



3.4 Kronecker Product Approximation 

In a time-dependent volatilities BS market the covariance matrix Emtv has time- 
dependent blocks. The multi-dimensional BM is the unique source of risk in the BS 
market and the generation of the trajectories of the 1-dimensional BM does depend 
on the volatilities. Based on these considerations, we propose to find a constant 
covariance matrix among the assets H, in order to approximate, in an appropriate 
sense, the matrix Emtv as a Kronecker product of R and H. In the following we 
illustrate the proposed procedure that we label Kronecker Product Approximation 
(KPA). Pitsianis and Van Loan [12] have proved the following proposition. 

Proposition 4. Suppose G G R mx ™ an d Gi e r™ix™i w ah m — m 1 m 2 and 
n = n\U2- Consider the problem of finding G\ E M™ 1 xni that realizes the minimum 

min || G - Gi ® G 2 \\%, (27) 

G 2 6l m i x "i 

where || • \\ F denotes the Frobenius norm. For fixed h — 1, . . . , to 2 and I = 1, . . . , n 2 
denote 1Z(G)hi the mi x n\ matrix defined by the rows h,h + m 2 , h + 2m2, . . . , h + 
(mi — l)m,2 and the columns I, I + n 2l I + 2fi2, . . . , I + {n\ — l)n2 of the original 
matrix G. The elements of G\ which gives {2ty are 

, s Tr (niG^Gi) , 
(G* 2 ) hl = ^(dGh h = l,...,m 2 ,l = l,...,n 2 , (28) 

where Tr indicates the trace of a matrix. 

In our setting, we have G = £mat, Gi = R and G2 = H . We note that for any 
i, j = 1, . . . , N, lZ(T,MN))ij is a N x N boomerang matrix. Moreover, given two 
general N x N boomerang matrices A and B, by direct computations we can prove 

N 

Tt(A t B) = Tr(AB) - £ (2(N - j) + 1) a n b n . (29) 
i=i 

Then we perform the PCA decomposition of R (8> H relying on the properties of 
the Kronecker product. However, if we use the PCA decomposition of the matrix 
F = R £§> H we do not get the required path. In order to produce the required 
trajectory we take 

Z - C KPA e - Cx MN (C^EnA^e (30) 

where Cs MJV and Cp are the CH matrices associated to Ejjjv and F, respectively, 

and EhA 1 ^ 2 is the PCA decomposition of F. The matrix C KPA produces the 

1/2 

correct covariance matrix; indeed, denoting P — E^A^ , we have 
C KPA (C KPA ) T = Cx MN (C F )-'PP T [{C F r'] T Cl MN = C^ MN Cl MN = E MJV 
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Table 1: Inputs Parameters 



Si(0) 




100, V* = 1 . . . , N 


K 


c 


{90,100,110} 


r 




4% 


T 




1 


<Ti(0) 




10% + ^±40% i = l...,N 


<7j( + 00) 




9% Vi = 1 . . . , N 


Ti 




1.5 W = l...,iV 


Pij 


c 


{0,40} i,j = l...,JV 



Table 2: Effective Dimensions. Time-dependent Volatilities 

p = 0% 

Ch PC A LT KPA 

dr > 1900 d T = 14 d T = 10 d T = 19 



p = 40% 

Ch PC A LT KPA 

d T > 1900 d T = 9 d T = 8 dr = 11 



because PP T = CpCp = i* 1 . Our fundamental assumption is that the principal 
components of Z are not so different from those of the normal random vector Z' 
whose covariance matrix is F. We expect that the KPA decomposition would 
produce an effective dimension higher than the effective dimension obtained by the 
PCA decomposition, but with an advantage from the computational point of view. 
Due to properties of the Kronecker product, Equation (|30)l becomes 

Z = C^ MN (C^ g C H l ) EnA^e, (31) 

where Cr and Ch are the CH matrices of R and H, respectively. This matrix mul- 
tiplication can be carried out quickly by block-matrices multiplication and knowing 
that, due to the Propositions [T] and [H C^ 1 is a sparse bi-diagonal matrix. 

4 Computing the Option Price 

We will now estimate the fair price of an Asian option on a basket of M = 10 un- 
derlying assets with N = 250 sampled points in the BS model with time-dependent 
volatilities having the following expression 

a,(t) =<Xi(0) exp(-i/Vi)+<7i(+oo), i=l,...M. (32) 

The parameters chosen for the simulation are listed in Table[T](<7i(0) is then (7,(0)- 
<7j(+oo)). We implement the numerical investigation in two parts: first we test 
the effectiveness of the path-generation constructions on dimension reduction and 
compute their computational times, and then we compare the accuracy of the 
simulation. 
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Table 3: Computational Times in Seconds 



Constant Volatilities 




p = 40% 



Ch PC A LT 
0.60 25.77 71.14 



Ch PC A LT 
0.59 25.55 71.02 



Time-dependent Volatilities 




p = 40% 



Ch PC A LT KPA 
0.62 565.77 71.65 28.25 



Ch PC A LT KPA 
0.62 568.55 71.20 28.33 



The Table [5] shows the effective dimensions obtained by all the path-generation 
methods (p = 0.99). The LT construction provides the lowest effective dimension, 
while the PCA decomposition performs almost as well as the LT approach for 
the correlation case only, and the KPA returns a slightly higher effective dimen- 
sion. The CH decomposition collects 98.58% and 98.70% of the total variance for 
dx ~ 2000 for the uncorrelated and correlated cases, respectively. To have a more 
accurate comparison, Table [3] displays the elapsed times computed in Sabino [8] 
using an ad hoc incremental QR algorithm for the LT and assuming constant 
volatilities each equal to o~i (0) of Table [TJ The computation is implemented in 
MATLAB running on a laptop with an Intel Pentium M, processor 1.60 GHz and 
1 GB of RAM. We compute 50 optimal columns for the LT technique. The CH al- 
gorithm for block boomerang matrices has almost the same cost as the one relying 
on the properties of the Kronecker product. As a consequence, the LT also requires 
almost the same computational cost, while the PCA needs a time almost 20 times 
higher because we can no longer rely on the properties of the Kronecker product. 
In contrast, the KPA has almost the same computational time as the PCA in the 
constant volatility case and is the best performing path-generation technique from 
the computational time point of view. We have applied Proposition [2] to imple- 
ment the PCA and computed the eigenvalues and eigenvectors of Emjv relying 
only on the sparse function of MATLAB. It is noteworthy to say that there ex- 
ist algorithms tailored for the computation of the eigenvalues and eigenvectors of 
tri-diagonal symmetric block matrices that can further reduce the computational 
time. 

In the second part of our investigation we launch a simulation in order to es- 
timate the Asian option price using 10 replications each of 8192 random points 
following the strategy in Imai and Tan [7]. We use different random generators: 
standard MC, LHS and RQMC generators. Concerning the computational times 
of the price estimation, the CPU ratio between LHS and RQMC is almost 1 while 
standard MC is 1.33 faster. Moreover, the LT construction needs a time that is 
almost 1/30 of the total computational time of the LHS or RQMC simulation. As 
a RQMC generator we use a Matousek scrambled version (see Matousek [T3]) of 
the 50-dimensional Sobol sequence satisfying Sobol's property A (see Sobol [To]). 
We pad the remaining random components out with LHS. This hybrid strategy is 
intended to investigate the effective improvement of the decomposition methods 
when coupled with (R)QMC. Indeed, it can be proven that LHS gives good vari- 
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Table 4: Estimated At-the Money Prices and Errors. 





p = 




P = 


40% 






Price 


RMSE 


Price 


RMSE 




Ch 


3.18200 


0.01300 


5.18900 


0.02600 


MC 


KPA 


3.12400 


0.01300 


5.19400 


0.02600 




PCA 


3.10600 


0.01300 


5.20100 


0.02600 




LT 


3.11100 


0.01300 


5.24200 


0.02600 






Price 


RMSE 


Price 


RMSE 




Ch 


3.12200 


0.00750 


5.20000 


0.01200 


LHS 


KPA 


3.12440 


0.00550 


5.20950 


0.00340 




PCA 


3.12010 


0.00540 


5.20070 


0.00320 




LT 


3.12200 


0.00290 


5.20090 


0.00120 






Price 


RMSE 


Price 


RMSE 




Ch 


3.11240 


0.00550 


5.19500 


0.01500 


RQMC 


KPA 


3.12240 


0.00086 


5.20080 


0.00054 




PCA 


3.12140 


0.00078 


5.20080 


0.00069 




LT 


3.12230 


0.00021 


5.20080 


0.00019 



ance reductions when the target function is the sum of one-dimensional functions 
(see Stein [H]). On the other hand, the LT method is conceived to capture the 
lower effective dimension in the truncation sense for linear combinations. As a 
consequence, we should already observe a high accuracy when running the simu- 
lation using LHS combined with LT. We expect the KPA technique to produce a 
suboptimal decomposition in the sense of ANOVA, with the advantage of a lower 
computational effort. Our setting is thought to test how large is the improvement 
given by all the factorizations. Tables 2] and [5] present the results of our investiga- 
tion. The prices in Table [4] are all in statistical agreement. Those obtained with 
the CH decomposition are almost not sensitive to the random number generator. 
KPA, PCA and LT all provide good improvements both for the LHS and RQMC 
implementations for all the strike prices. The LT has an evident advantage com- 
pared to the PCA and KPA constructions in the uncorrelated case. In contrast, we 
observe that the KPA and PCA-based simulations give almost the same accuracy, 
both assuming uncorrelated and correlated asset returns. Considering the total 
computational cost and accuracy we observe that the KPA performs better than 
the standard PCA. Moreover, all these constructions can be employed in stochas- 
tic and local volatility models that are based on the mixture of multi-dimensional 
dynamics for basket options as done in Brigo et al. [3J. 



5 Computing the Sensitivities 

In the financial jargon, a Greek is the derivative of an option price with respect 
to a parameter. A Greek is therefore a measure of the sensitivity of the price 
with respect to one of its parameters. The deltas (A's) are the components of 
the gradient of the discounted expected outcome of the option with respect to 
the initial values of the assets. The problem of computing the Greeks in finance 
has been studied by several authors. In the following we extend the methodology 
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employed by Montero and Kohatsu-Higa pQ , based on the use of Malliavin Calculus, 
to the multi-assets case. The main difficulties of this extension lie in the fact 
that the assets are now correlated and the formulas in Montero and Kohatsu- 
Higa [1 cannot be directly extended to the multi-dimensional case. Moreover, the 
localization technique, introduced by Fournie et al [9], should generally control 
all the components of the multi-dimensional BM to improve the accuracy of the 
estimation. We write the dynamics (Q} with respect to a M-dimensional BM W(t) 
with uncorrelated components 



M 



dSi(t) = rSi(t)dt + Si{t)ai{t) V (t)dW m (t) * 



m— 1 



,Af, 



(33) 



where Y,m=i a im^km = Pik and we have defined <7 im (t) 

The Malliavin calculus is a theory of variational stochastic calculus and provides 
the mechanics to compute derivatives and integration by parts of random variables 
(see Nualart [T7] for more on Malliavin Calculus). 

Denote by D\ , . . . , Df 1 the Malliavin derivatives with respect to the compo- 
nents of W(t), while S sk = X)m=i ^nt represents the Skorohod integral with 5^ 
indicating the Skorohod integral on the single W m (t). The domains of the Malliavin 
derivatives and the Skorohod integral are denoted by D 1 ' 2 and dom((5 sk ), respec- 
tively, while 5 Kr indicates the Kronecker delta. We prove the following proposition. 

Proposition 5. Denote x = S(0), and G k the partial derivative 



dm(T) _ Ef=iWkjS k (tj) 



G k = 



dxk 



Xk 



k = 1, . 



(34) 



where m(T) — Y^,jLi w kjSk(tj). Knowing that a(T) S B 1,2 , the k-th delta (the k-th 
component of the gradient) is 



dx k 



e- rT E [a'(T)G k ]=e- rT E 



M 



a(T) ^ 6: 



Sk 



(G k u m ) 



(35) 



where u = (u\, 
dom(<5 sk ) and 



,u M ) € dom(<5 bk ) 



0*1 1 



%m (^) 

EH=ifo^(s)D^m(T)ds 

M 



G dom(<5 Sk ), G fe u e 



[ z h {s)D h s 7n(T)ds / 0, a.s. 



Proof. Compute 



D h s a{T) = a' \T)D h s m(T) h = 1, 



,M. 



Multiply the above equation by Zhit) - so that z 6 dom(i5 Sk ) 
sum for all h = 1, . . . , M and integrate: 



(36) 

and by G k ; then 



M 

E 

h=l 



G k z h (s)D h s a(T)ds 



M 

E 

h=l 



G k z h {s)a\T)D h s m{T)ds. (37) 
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Due to the definition of u and to the fact that a'(T)Gk does not depend on s we 
can write 



M T 

a\T)G k = Y, / u m (s)G k D™a(T))ds. k = 1, 

m=l J 

Finally compute the expected value of both sides of (|38|) 

M T 



,M 



(38) 



V / u m (s)G k D™a(T)ds 



E[a'(T)G fc ] =E 
so that by duality 

A k = E [a(T))S sk (G k u)] k = 1, . 
and this concludes the proof. 



(39) 

(40) 
□ 



Proposition [5] allows a certain flexibility in choosing either the process u, or 
better z. We consider z% — a k S^; h,k = 1,...,M, a k — l,Vfc. Namely in 
order to compute the fc-th delta we consider only the fc-th term of the Skorohod 
integral reducing the computational cost. In particular, this choice is motivated 
by the fact that in this way the localization technique needs to control only <S^ k (-) 
and then only the fc-th component of W(i). Then we define L k and calculate for 
k = l,...,M 

M N 



D k G k ds 



D k L k ds 



f D k s m(T)ds = y"f>ySifo) [\ ik (s)ds, (41) 
J o i= i j= i Jo 

2_, w ]kS k {t :j ) I a kk (s)ds = y] w jkSk(tj) a k (s)ds, (42) 
<r,,S,:l ,) ^ <r ik ds 



(43) 



and hence 



A fc =E 



,M. 



(44) 



Due to the properties of the Skorohod integral we have for k = 1, . 



S k 



Gfe 



Gfc 
L K 



W k {T) - J? | u, 



L k / ^'Gfcds - G k / ^Lfcds , (45) 



With another choice of z, for instance Zh — olhi A k would depend linearly on the 
whole Af-dimensional BM, making the localization technique less efficient. 

We investigate the applicability of the RQMC approach to estimate the ex- 
pected value in Equation (|44|) for fc = 1, ...,M. We assume the same input 



parameters as in Section [4] and generate the trajectories (the values Si(tj),i — 
1, . . . , M, j — 1, . . . , N) in Equations (|41M3| as in Section |H Moreover, we con- 
sider ai m as the elements of the CH matrix associated to pi m , i,m — l,...,M. 
Table [5] compares the deltas obtained with RQMC only with the same number of 
scenarios as in Section [4] We apply the same LT construction used to estimate the 
price of the option and not the one for the integrand function in Equation 
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Table 5: At-the-money estimated A's (1(T 2 ) and errors (10~ 4 ) with RQMC. 



p = 0% 

LT KPA PCA CH 



A 


RMSE 


A 


RMSE 


A 


RMSE 


A 


RMSE 


6.1832 


0.80 


6.1820 


1.10 


6.1808 


0.86 


6.2060 


1. 


.50 


6.2024 


0.75 


6.2126 


0.90 


6.2016 


0.86 


6.2250 


1 


.10 


6.2305 


0.85 


6.2340 


0.87 


6.2341 


0.99 


6.2530 


1. 


,20 


6.2667 


0.75 


6.2701 


0.82 


6.2699 


0.91 


6.2830 


1. 


,40 


6.3081 


0.60 


6.3133 


0.92 


6.3093 


0.96 


6.3270 


1. 


.20 


6.3569 


0.55 


6.3595 


0.97 


6.3598 


0.83 


6.3750 


1. 


.10 


6.4107 


0.50 


6.4141 


0.93 


6.4103 


0.78 


6.4329 


1. 


,20 


6.4709 


0.50 


6.4744 


0.91 


6.4677 


0.84 


6.4920 


1. 


.40 


6.5338 


0.50 


6.5390 


0.93 


6.5325 


0.93 


6.5530 


1. 


,30 


6.6001 


0.65 


6.6060 


0.78 


6.6000 


0.96 


6.6120 


1. 


.10 



p = 40% 

LT KPA PCA CH 



A 


RMSE 


A 


RMSE 


A 


RMSE 


A 


RMSE 


5.47830 


0.056 


5.48410 


0.110 


5.48050 


0. 


130 


5.46767 


1.100 


5.53510 


0.062 


5.54050 


0.110 


5.53740 


0. 


140 


5.52457 


1.200 


5.59430 


0.054 


5.60020 


0.110 


5.59660 


0. 


130 


5.58730 


1.200 


5.65440 


0.062 


5.66120 


0.120 


5.65680 


0. 


130 


5.64031 


1.000 


5.71680 


0.075 


5.72330 


0.130 


5.71790 


0. 


140 


5.70994 


1.100 


5.78130 


0.082 


5.78850 


0.120 


5.78410 


0. 


120 


5.77038 


1.300 


5.84840 


0.077 


5.85320 


0.096 


5.85060 


0. 


120 


5.83233 


1.200 


5.91560 


0.082 


5.92110 


0.110 


5.91790 


0. 


130 


5.90019 


1.100 


5.98490 


0.052 


5.99070 


0.093 


5.98680 


0. 


120 


5.97059 


1.000 


6.05470 


0.059 


6.06050 


0.110 


6.05680 


0. 


130 


6.04613 


1.200 
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This would not seem to be the optimal choice, but if we would have applied the LT 
for the integrand function in Equation (|4"4")l M = 10 decomposition matrices (one 
for each delta) would be required. This setting would have increased the CPU 
time to obtain the LT to at least 1/3 (even higher due to the larger number of 
terms to compute) of the total time making the estimation less convenient. Table 
[5] shows that the PCA, LT and KPA approaches perform almost equally in terms 
of RMSEs, with the LT giving only slightly better results in the uncorrelated case. 
In terms of computational cost the KPA performs better than the PCA. The CH 
construction displays RMSEs that are even 10 times higher. As explained before, 
g can be considered a good approximation for the payoff function in Equation (|4~4")l 
but in the Malliavin expression a(T) is multiplied by a random weight that depends 
on the Gaussian vector Z. In contrast, the PCA and the KPA concentrate most 
of the variation in the first dimensions of Z. This is the explanation of the almost 
equal accuracy of the LT, the PCA and the KPA. 

6 Conclusions 

We have considered the problem of computing the fair price and the deltas of high- 
dimensional Asian basket options in a BS market with time-dependent volatilities. 
In order to extend the QMC superiority to high dimensions it is necessary to 
employ path-generation techniques with the main purpose to reduce the nominal 
dimension. The LT and the PCA constructions try to accomplish this task by the 
concept of ANOVA. In the case of time-dependent volatilities in the BS economy 
the computational cost of the LT and the PCA cannot be reduced relying on the 
properties of the Kronecker product and the computation is more complex. We 
have presented a new and fast CH algorithm for block matrices that remarkably 
reduces the computational burden making the LT construction even more conve- 
nient than the PCA. We have introduced a new path-generation technique, named 
KPA, that in the applied setting, is as accurate as the PCA and is even more con- 
venient with respect to the computational cost. In addition, we proved that the 
KPA enhances RQMC for the estimation of the fair price and the calculation of the 
deltas of Asian basket options in a BS model with time-dependent volatilities. In 
this setting the KPA provides the same accuracy of the LT in the case of correlated 
asset returns and in the estimation of the deltas. Finally, concerning the computa- 
tion of the sensitivities, we have extended the procedure adopted by Montero and 
Kohatsu-Higa [T], based on the Malliavin Calculus, to the multi-assets case. All 
these results can be easily applied to stochastic and local volatility models that are 
based on the mixture of multi-dimensional dynamics for basket options, as done in 
Brigo et al. [2]. 
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